EXAMPLE #1: Closed Economy Sidrauski Model. See paper for details.
Manoj Atolia: September 18, 2007.
Works with both Mathematica 6.0 and 5.2.

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[FindMinimum::lstol];Off[FindMinimum::reged];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* This cell sets the precision of calculations. *)
{precision = 100,$MinPrecision=digits=precision};

Initialization

(* Set the directory where Mathematica will store intermediate results. *)
SetDirectory["c:"];

(* ALWAYS ENTER NUMERICAL VALUES AS FRACTIONS AND NOT AS DECIMALS TO ENSURE REQUISITE PRECISION IN CALCULATIONS *)   
{τ=(1/4),ρ=(1/10),β=(1/2),δ=(1/20)}//N                                   (* Parameters of the model *)

{(Subscript[bbyY, o])=(25/100),(Subscript[mbyY, o])=(15/200),(Subscript[π, o])=(1/2),(Subscript[θ, o])=(1/2)}//N        (* Data for calibration *)

{(Subscript[z, o])=1,(Subscript[k, o])=1}//N                                                                          (* Normalization *)

{(Subscript[z, f])=(6/5)(Subscript[z, o]),ω=(1/5)}//N                                                                 (* Shock - Exogenous rise in productivity *)

Steady State

(* Define utility function and production function *)
U[c_,m_]:=((1 )/(1-(1/τ)))((c^(1-(1/τ)))+g (m^(1-(1/τ))));Y[k_,z_]:=z((((Subscript[a, 0])+(Subscript[a, 1])(k^((β-1)/β))))^(β/(β-1)));(Subscript[Y, 1])[k_,z_]:=Evaluate[D[Y[k,z],k]];
(* Avoid using symbols for subscripts for representing derivatives. Thus, (Subscript[Y, 1])[k,z] is (Subscript[Y, K])[k,z] of the paper.*)
(* Since we are using symbols o and f as subscripts for initial and final steady states avoid using o and f as names of variables or function. *)

(* Calibrate and compute the original steady state *)
(Subscript[Y, o])=Y[(Subscript[k, o]),(Subscript[z, o])];{(Subscript[c, o])=(Subscript[Y, o])-δ (Subscript[k, o]),(Subscript[m, o])= (Subscript[mbyY, o]) (Subscript[Y, o]),b=(Subscript[bbyY, o]) (Subscript[Y, o]),L=(Subscript[π, o])(Subscript[m, o])-ρ b};g=((Subscript[π, o])+ρ)(((((Subscript[m, o])/(Subscript[c, o]))))^(1/τ));{(Subscript[a, 0]),(Subscript[a, 1])}=.;{(Subscript[a, 0]),(Subscript[a, 1])}={(Subscript[a, 0]),(Subscript[a, 1])}/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, o]),(Subscript[z, o])] == ρ+δ,(Subscript[θ, o])==(((ρ+δ)(Subscript[k, o]))/(Subscript[Y, o]))},{(Subscript[a, 0]),1},{(Subscript[a, 1]),1},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];

(* Compute the new steady state *)
(Subscript[k, f])=.;(Subscript[k, f])=(Subscript[k, f])/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, f]),(Subscript[z, f])]==ρ+δ},{(Subscript[k, f]),(Subscript[k, o])},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];(Subscript[Y, f])=Y[(Subscript[k, f]),(Subscript[z, f])];(Subscript[c, f])=(Subscript[Y, f])-δ (Subscript[k, f]);
{(Subscript[m, f]),(Subscript[π, f])}=.;{(Subscript[m, f]),(Subscript[π, f])}={(Subscript[m, f]),(Subscript[π, f])}/.FindRoot[{g(((((Subscript[m, f])/(Subscript[c, f]))))^(-1/τ))==(Subscript[π, f])+ρ,L==(Subscript[π, f])(Subscript[m, f])-ρ b},{(Subscript[m, f]),(Subscript[m, o])},{(Subscript[π, f]),(Subscript[π, o])},AccuracyGoal->digits,WorkingPrecision->digits];
{{(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[Y, o]),b,L,(Subscript[z, o]),(Subscript[π, o]),(Subscript[c, o])/(Subscript[m, o]),(Subscript[m, o])/(Subscript[Y, o])},{(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f]),(Subscript[Y, f]),b,L,(Subscript[z, f]),(Subscript[π, f]),(Subscript[c, f])/(Subscript[m, f]),(Subscript[m, f])/(Subscript[Y, f])}}//N//MatrixForm

Global Nonlinear Saddle Path

(* Specify the nonlinear system with state variables at the end *)
{mdot = L+r b - pi m; 
cdot = τ c (r-ρ); 
kdot = Y[k,z]-c-δ k; 
zdot = ω ((Subscript[z, f])-z);
system  ={mdot,cdot,kdot,zdot};};

(* Solve for the endogenous variables as functions of the variables in the system. For consistency maintain same order of variables in arguments for functions of endogenous variables as in the system *)
{r[k_,z_]:=(Subscript[Y, 1])[k,z]-δ;
pi[m_,c_, k_,z_]:=g((((m/c)))^(-(1/τ)))-r[k,z];};

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
{  tol=(10^(-5));                                                (* minimum distance of acceptable trajectory from initial steady state *);
  precision = MachinePrecision;                                                              (* starting value of precision of calculations *);
  epsilon= (10^(-2))Norm[{(Subscript[k, f])-(Subscript[k, o]),(Subscript[z, f])-(Subscript[z, o])}];               (* radius of the ϵ-ring*);};

?solve

(* MAIN CELL - Solving for the global nonlinear saddlepath. *)
Timing[{gtmin,gminsol}=solve[system,{t},{m,c,k,z},{(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[z, o])},{(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f]),(Subscript[z, f])},{{r,{k,z}},{pi,{m,c,k,z}}},{tol,precision,epsilon}]]

(* Plot the final solutions. *)
tplot = 65;
mPlot=Plot[{If[t<0,(Subscript[m, o]),m[gtmin-t]]/.gminsol,(Subscript[m, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",m},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
cPlot=Plot[{If[t<0,(Subscript[c, o]),c[gtmin-t]]/.gminsol,(Subscript[c, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",c},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
kPlot=Plot[{If[t<0,(Subscript[k, o]),k[gtmin-t]]/.gminsol,(Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",k},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
zPlot=Plot[{If[t<0,(Subscript[z, o]),z[gtmin-t]]/.gminsol,(Subscript[z, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",z},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
YPlot=Plot[{If[t<0,(Subscript[Y, o]),Y[k[gtmin-t],z[gtmin-t]]]/.gminsol,(Subscript[Y, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",Y},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
IPlot=Plot[{If[t<0,δ (Subscript[k, o]),Y[k[gtmin-t],z[gtmin-t]]-c[gtmin-t]]/.gminsol,δ (Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"","I"},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
rPlot=Plot[{If[t<0,ρ,r[k[gtmin-t],z[gtmin-t]]]/.gminsol,ρ},{t,-tplot/10,tplot},PlotRange->All,AxesOrigin->{-tplot/10,ρ},AxesLabel->{"",r},DisplayFunction->Identity];
piPlot=Plot[{If[t<0,(Subscript[π, o]),pi[m[gtmin-t],c[gtmin-t],k[gtmin-t],z[gtmin-t]]]/.gminsol,(Subscript[π, f])},{t,-tplot/10,tplot},PlotRange->All,AxesOrigin->{-tplot/10, (Subscript[π, f])},AxesLabel->{"",π},DisplayFunction->Identity];
statePlot=ParametricPlot[{k[gtmin-t],z[gtmin-t]}/.gminsol,{t,0,tplot},PlotRange->All,AxesLabel->{k,z},DisplayFunction->Identity];

If[$VersionNumber==6.0,
Show[GraphicsGrid[{{mPlot,cPlot,kPlot},{zPlot,YPlot,IPlot},{rPlot,piPlot,statePlot}}]],
Show[GraphicsArray[{{mPlot,cPlot,kPlot},{zPlot,YPlot,IPlot},{rPlot,piPlot,statePlot}}]]]

(* Tabulate the final solutions. *)
Labels={"m","c","k","z","Y","I","r","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[z, o]),(Subscript[Y, o]),δ (Subscript[k, o]),ρ,(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f]),(Subscript[z, f]),(Subscript[Y, f]),δ (Subscript[k, f]),ρ,(Subscript[π, f])}//N;
(Dat=Table[({tt,m[t],c[t],k[t],z[t],Y[k[t],z[t]],Y[k[t],z[t]]-c[t],r[k[t],z[t]],pi[m[t],c[t],k[t],z[t]]}/.{t-> gtmin-t}/.tt->t/.gminsol)[[1]],{t,0,65,1/4}])//MatrixForm;
(Dat=Join[{Join[{"Initial Steady State"},InitialSteadyState],Join[{"Final Steady State"},FinalSteadyState],Join[{Time},Labels]},Dat])//N//MatrixForm

The solve 2D Fast Package

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[FindMinimum::lstol];Off[FindMinimum::reged];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* Interface *)
BeginPackage["solve`"];

EndPackage[];

(* Usage messages for the exported functions and the context itself *)
solve::usage := 
"solve[sys, {t}, vars, varso, varsf, varsendo, init] solves a system of first-order differential equations. \n
sys   = {xdot, ydot,...} lists the right hand side of the system of first-order differential equations (with state variables at end).\n
vars  = {x, y, ...} is the list of variables in the system in the same order as the differential equations in sys.\n
varso = {(Subscript[x, o]), (Subscript[y, o]), ...} is the list of original steady values of variables in the system in the same order as vars.\n
varsf = {(Subscript[x, f]), (Subscript[y, f]), ...} is the list of final steady values of variables in the system in the same order as vars.\n
varsendo = {{w,{x,y, ...}}, ...} is the list of the endogenous variables {w, ...} where w is function of {x,y, ...} where {x,y,...} are in the order in which they appear as arguments of w. If no endogenous variables enter {{,{}}}.\n
init   = {tol, precision, epsilon} sets the required accuracy to tol, starting precision to precision, and radius of initial ϵ-ring to epsilon."

(* Error messages for the exported objects *)
{solve::sysDESys = "You called `1` with first argument `2` (the system of differential equations) that has variables other than those defined in `3` and `4`.";
solve::sysLength = "You called `1` with lengths of `2`, `3`, `4` and `5` not equal.";
solve::vars = "You called `1` with third argument `2` that is not a list of variables.";
solve::origVals = "You called `1` with fourth argument `2` that is not a list of numbers.";
solve::finVals = "You called `1` with fifth argument `2` that is not a list of numbers.";
solve::endoVars = "You called `1` with sixth argument `2` that is not a list of lists of the form {w, {x,y, ...}} where w is an endogenous variable and {x,y, ...} are the variables in the system as defined in `3`.";
solve::initVals = "You called `1` with last argument `2` that is not a list of numbers.";};

(* Implementation *)
Begin["solve`Private`"];

(* Definition of the exported functions *)
solve[asys_List,{at_},avars_List,avarso_List,avarsf_List,avarsendo_List,ainit_List]:=Block[{sys,t,vars,varso,varsf,varsendo,init,tol,precision,epsilon,Tmax,firstring,
firstrun,box,system,j,ntraj,firstit,ini,pmin,sol,tsol,dis,tdis,gmindis,jmin,j1,j2,doneflag,T,gtmin,gminsol},

{sys,{t},vars,varso,varsf,varsendo,init}={asys,{at},avars,avarso,avarsf,avarsendo,ainit};{tol,precision,epsilon}=init;
Tmax=4000;firstring = 01;firstrun = 01;

(* Check input *)
checkInput[];

(* Make the system from input data and create box over which to solve the system *)
system=makeSys[];
box:=Apply[Or,Join[MapThread[First[Map[(#1),{t}]]<((Min[#2,#3])/100)&,{vars,varso,varsf}],MapThread[First[Map[(#1),{t}]]>5 Max[#2,#3]&,{vars,varso,varsf}]]];

(* SOLVE FOR THE GLOBAL NONLINEAR SADDLEPATH *)
	While[True,
	Label[precbegin];
	$MinPrecision=digits=precision;ntraj=20;j1=0;j2=ntraj;gmindis=(10^7);firstit=1;

	Label[begin];
	For[j=j1,j≤ j2,j++,
		ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==(varsf+Join[Table[0,{i,1,Length[vars]-2}],{epsilon Cos[2π j/ntraj],epsilon Sin[2 π j/ntraj]}]),List];

		If[((IntegerPart[j/2]-j/2)≠0||firstit==1),
			If[firstit==0&&Max[Abs[{ini[[-2]]/pmin[[1]]-1,ini[[-1]]/pmin[[2]]-1}]]<(10^(-precision)),
				precision=precision+4;
				Print["Suspecting loss in precision."];Print["Increased 'precision' by 4 digits to ",precision," digits."];Goto[precbegin]];

			sol=NDSolve[SetPrecision[{system,ini},precision],vars,{t,0,Tmax},
				Method->{"EventLocator", "Event"->SetPrecision[box,precision],"EventAction":> Throw[Null, "StopIntegration"],"EventLocationMethod"->"StepBegin"},
					MaxSteps->100000,WorkingPrecision->digits];tsol=sol[[1,1,2,1,1,2]];{dis,tdis}=distance[sol,tsol];
	
			If[firstit==1,
				Print["Trajectory = ",j,"     Duration of Solution = ",tsol//N,"     Minimum Distance = ",dis,"     Time of Closest Approach = ",tdis]];

			If[dis≤gmindis,gmindis=dis;gtmin=tdis;jmin=j;gminsol=sol ;pmin={ini[[-2]],ini[[-1]]};
				If[gmindis<tol/3,Print["Target Tolerance = ",tol/3//N,"   Tolerance Achieved = ",gmindis];  Goto[end]];];];
	];

	firstit=0;
	printAndPlot[];
	{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};Goto[begin];

	Label[end];
	If[firstring==1 && firstrun==1, T=IntegerPart[2 gtmin/3]>>T,T=<<T];

	(* Error Checking *)
	{firstring,firstrun, epsilon, precision, doneflag}=errorCheck[];
	If[doneflag==1,Break[]];
];

(* Check the duration for which the solution is accurate *)
durationCheck[];

{gtmin,gminsol}];

checkInput[___]:=Block[{},

If[Not[Length[sys]==Length[vars]==Length[varso]==Length[varsf]],Message[solve::sysLength,"'solve'","'sys'","'vars'","'varso'","'varsf'"];Abort[]];

If[Not[Length[vars]==Length[Cases[vars,_Symbol]]],Message[solve::vars,"'solve'","'vars'"];Abort[]];

If[Not[Length[varso]==Length[Cases[varso,_?NumericQ]]],Message[solve::origVals,"'solve'","'varso'"];Abort[]];

If[Not[Length[varsf]==Length[Cases[varsf,_?NumericQ]]],Message[solve::finVals,"'solve'","'varsf'"];Abort[]];

If[Not[Length[varsendo]==Length[Cases[varsendo,{_Symbol,_List}]]]||Not[Length[Flatten[varsendo]]==Length[Cases[Flatten[varsendo],_Symbol]]],Message[solve::endoVars,"'solve'","'varsendo'","'vars'"];Abort[];];

If[Not[Length[init]==Length[Cases[init,_?NumericQ]]],Message[solve::initVals,"'solve'","'init'"];Abort[]];

];

makeSys[___]:=Block[{system},

system= -sys/. Thread[MapThread[(#1)&,Transpose[varsendo]]->
				MapThread[ToExpression[StringJoin[ToString[ #1],"[",ToString[#2],"]"]] /.List->Sequence&, Transpose[varsendo]],
					List];

If[Length[Cases[system/.Thread[vars->varsf,List],_?NumberQ]]≠Length[vars],
	Message[solve::sysDESys,"'solve'","'sys'","'vars'","'varsendo'"];Abort[];];

system=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==
			(system=system/.Thread[vars->Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars],List]),
				List];

system];

printAndPlot[___]:=Block[{},

Print["Extreme Left Trajectory = ",j1,"   Closest Trajectory = ",jmin,"   Extreme Right Trajectory = ",j2,"    No. of Subdivisions of the ϵ-ring = ",ntraj];

Print["Current Minimum Distance = ",gmindis,"  Time of Closest Approach = ",gtmin];

Print[ParametricPlot[{{vars[[-2]][t],vars[[-1]][t]}/.gminsol,{varso[[-2]]+(1/100) (varsf[[-2]]-varso[[-2]]) Sin[t],varso[[-1]]+(1/100) (varsf[[-1]]-varso[[-1]]) Cos[t]}},
					{t,0,gtmin},PlotRange->All]];

];

distance[sol_,tsol_]:=If[firstit==1,
Block[{dis,dislast,steps,tbeg,tend,intsz,traj,tdis},

dis= (10^7);dislast=(10^8);steps=Max[10Floor[tsol],450] ;tbeg=0;tend=tsol;intsz=tend-tbeg;

While[ (dislast-dis)>tol (10^(-5)),

dislast=dis;
traj=Table[Norm[Flatten[{vars[[-2]][t]-varso[[-2]],vars[[-1]][t]-varso[[-1]]}/.sol],2],{t,tbeg,tend,(intsz/steps)}];
dis=Min[traj];tdis=tbeg+(Position[traj,dis][[1,1]]-1)(intsz/steps);
tbeg=Max[tdis-(intsz/steps),0];tend=Min[tdis+(intsz/steps),tsol];intsz=tend-tbeg;steps=40];

{dis,tdis}],

Block[{dis,tdis,tbeg,tend,min},

tbeg=0;tend=tsol;

min=FindMinimum[First[Norm[Flatten[{vars[[-2]][t]-varso[[-2]],vars[[-1]][t]-varso[[-1]]}]]/.sol],{t,9tsol/10,tbeg,tsol},WorkingPrecision->digits,AccuracyGoal->5-Log[10,tol]];

{dis,tdis}={min[[1]],min[[2,1,2]]}]];

errorCheck[___]:=Block[{diffprec,diffeps, satisprec, notsatisprec, notsatiseps, increaseprec,decreaseps, step2, step3, congratulations, done,gtmin1,gminsol1},

doneflag=0;

(* Messages *)
diffprec:=Print["Difference in values of variables due to latest increase in precision = ", sysDiff[0,T]//N];
diffeps:=Print["Difference in values of variables due to latest decrease in epsilon = ",sysDiff[0,T]//N];

satisprec= "THE PRECISION IS SATISFACTORY. The variables changed by less than 'tol' due to latest increase in 'precision' for current ϵ-ring.";
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than 'tol' due to latest increase in 'precision' for current ϵ-ring.";
increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];

decreaseps:=Print["Decreased 'epsilon' by a half to ",N[epsilon,precision],"."];
notsatiseps="The epsilon rings are NOT tight enough. The variables changed by greater than 'tol' due to latest decrease in 'epsilon'.";

step2="Starting Step 2: Checking for the precision of the candidate solution.";
step3="Starting Step 3: Checking for the accuracy of the candidate solution.";

congratulations= "CONGRATULATIONS! The variables changed by less than 'tol' due to latest decrease in 'epsilon'.";
done = "You are done! Hope to see you in another model. Continue with the next cell.";

Which[  (* possible to get rid of step 2 like the 3D case *)

firstring == 1 && firstrun == 1,
	{gminsol,gtmin}>>solution1;Print[step2];precision=precision+4;increaseprec;firstrun=0,firstring == 1 && firstrun == 0,
	{gminsol1,gtmin1}=<<solution1;{gminsol,gtmin}>>solution1;diffprec;
	If[Max[sysDiff[0,T]]>tol,
		Print[notsatisprec];precision=precision+4;increaseprec,
		Print[satisprec];Print[step3];firstring=0;firstrun=1;epsilon=epsilon/2;decreaseps],

firstring == 0  && firstrun == 1,
	{gminsol,gtmin}>>solution2;{gminsol1,gtmin1}=<<solution1;
	If[Max[sysDiff[0,T]]>tol,
		firstrun=0;precision=precision+4;increaseprec,
		diffeps;Print[congratulations];Print[done];doneflag=1 ],

firstring == 0  && firstrun == 0,
	{gminsol1,gtmin1}=<<solution2;{gminsol,gtmin}>>solution2;diffprec;
	If[Max[sysDiff[0,T]]>tol,
		Print[notsatisprec];precision=precision+4;increaseprec,
		Print[satisprec];{gminsol1,gtmin1}=<<solution1;diffeps;
		If[Max[sysDiff[0,T]]>tol,
			{gminsol,gtmin}>>solution1;Print[notsatiseps];firstrun=1;epsilon=epsilon/2;decreaseps,Print[congratulations];
			Print[done];doneflag=1]]
];

{firstring,firstrun,epsilon,precision,doneflag}];

sysDiff[t1_?NumberQ,t2_?NumberQ]:=Table[Max[Abs[Table[((vars[[i]][gtmin-t]/.gminsol)-(vars[[i]][gtmin1-t]/.gminsol1))[[1]],{t,t1,t2,1/10}]]],{i,1,Length[vars]}];

durationCheck[___]:=Block[{gminsol1,gtmin1,j},

{gminsol1,gtmin1}=(<<solution1)[[{1,2}]]; {gminsol,gtmin}=(<<solution2)[[{1,2}]]; T=<<T;

For[j=T+1,j≤IntegerPart[gtmin1],j++,
	If[Max[sysDiff[j-1,j]]>tol,
		Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years."];Break[]]];

If[j==IntegerPart[gtmin1]+1,Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years."]];Print["The solution upto ",T," years is accurate to within ",Max[sysDiff[0,T]]//N,". Continue with next cell."]

];

End[];

(* Epilog *)
Protect[Evaluate[$Context<>"*"]];